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Abstract 

We investigate numerically the behaviour of damage spreading in a Kauffman cellular automaton 
with quenched rules on a dynamical (f> 3 graph, which is equivalent to coupling the model to dis- 
cretized 2D gravity. The model is interesting from the cellular automaton point of view as it lies 
midway between a fully quenched automaton with fixed rules and fixed connectivity and a (soluble) 
fully annealed automaton with varying rules and varying connectivity. In addition, we simulate the 
automaton on a fixed </> 3 graph coming from a 2D gravity simulation as a means of exploring the 
graph geometry. 
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Cellular automata have found applications in many area of physics and other sciences. One particularly 
intriguing idea was put foward by Kauffman in the 1960s fil|, who suggested that a random Boolean 
network might serve as a model for cell differentiation. This predated the more recent explosion of 
papers on the subject that was largely inspired by the work of Wolfram H. In the Kauffman model each 
oi i = 1,2,. . . , N lattice sites contains a spin which can take the values ±1 (or equivalently a Boolean 
variable that can take the values 0, 1) and the spins evolve according to 

a i (t+l) = f i (a jl (t),<T j2 (t),...,a jh (t)) (1) 

where the functions fi are chosen independently and each site i has K specified inputs sites j± (i) , . . . ,jk(i)- 
The original biologically motivated model of Kauffman chose the fi randomly at t = with a probability 
p = 1/2 for fi = +1 and a probability 1 — p = 1/2 for /j = — 1. for each site and also chose the 
J2(*)i ■ ■ ■ j jk(i) inputs randomly from the N lattice sites at t = 0. This version of the model 
thus displays quenched randomness as both the functions fi and the input sites are fixed and has been 
called "infinite dimensional" because of the non-locality of the inputs. The model must display periodic 
behaviour for finite N as there are only 2 N possible spin configurations, so one can enquire about the 
length and number of limit cycles as well as quantities such as the Hamming distance between different 
configurations 1 and 2, Dn{t)i defined as 

1 N 

^w^EWw-^w) 2 ( 2 ) 

where a 1 is a spin in configuration 1 and a 2 is a spin in configuration 2. The application to genetics 
arises from identifying the stable limit cycles with different stable genotypes. 

It was found (numerically) (^] that for K < 2 the length of cycles was small, of order ^{^N 1 / 2 , 
whereas for K > 3 it increased as exp(a(iV)iV), which implied a critical K c between 2 and 3, with 
corresponding singular behaviour in a{K) and (3(K). The K < 2 phase was called a frozen phase and 
had 

lim D 12 (t) = (3) 

t— >oo 

for any starting configurations 1 and 2, whereas the K > 2 phase was chaotic, having 

lim D 12 (t) = d final (4) 

t — >oo 

with the final distance dfi na i independent of the initial distance. A simple annealed approximation to 
the model was solved by Derrida and Pomeau [[| and captured these features giving a critical value of 
K = 2. One can think of implementing the approximation in a simulation by changing both the functions 
and the inputs and each time step, which amounts to neglecting all the correlations. 

In || Derrida and Stauffer investigated a finite dimensional variant of the model in which the K 
inputs were chosen from the neighbours of a site on various lattices. They varied the probability p, (1 —p) 
for choosing fi = +1, (— 1), which can be shown to give the condition 

2Kp c (l-p c ) = l (5) 

for the transition between the frozen and chaotic phases, where p is now used as a control parameter rather 
than K. Derrida and Stauffer pointed out that the Kauffman model with varying fi on a d dimensional 
lattice was equivalent (as far as distances and overlaps between configurations were concerned) to d + 
1 dimensional directed percolation || on the appropriate lattices. The qualitative properties of this 
annealed model were similar to the infinite dimensional case, although the numerical values of thresholds 
and exponents were different. 

They also investigated numerically a version with fixed fi and found rather different behaviour from 
the infinite dimensional case, with the final distance in the chaotic phase depending on the initial distance 
(cf equ.(4)). The behaviour is summarized as: 

lim L>i 2 (oo) = Q(p) > 0, if p > p c 

D 12 (0)— >0 
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lim D 12 {oo) = 0, if p < p c 

1 F? D ^ 2 ^\ = X(P)> */ P<Pc (6) 
Di 2 (o)^o D 12 {0) 

where the susceptibility x(p) diverges at p = p c . The strategy adopted in || was to look at the time 
variation of the distance between two configurations, where one of the initial configurations contained 
some "damaged" spins that were different to the other reference configuration. It is possible to think 
of the transition from the frozen phase to the chaotic phase as damage percolating through the lattice, 
which allows another estimate of p c if one neglects spin-spin correlations, namely 

2p c (l -Pc) = x c (7) 

where x c is the bond percolation threshold for the lattice. Numerical simulations of the quenched model 
on various lattices 0, |8| showed that frozen and chaotic phases existed on square, triangular and various 
cubic lattices, whereas the honeycomb (hexagonal) lattice possessed only a frozen phase. It was found 
that estimates of p c were subject to strong finite size effects, with only the square lattice value being 
pinned down with precision. In addition, various interesting quantities such as the fractal dimension of 
the cluster of damaged spins at threshold were measured. The analogy with percolation for the quenched 
model seems to be rather less solid than the annealed version, as various thresholds and fractal dimensions 
differ from the percolative case |§ . 

The work of |], |(| |?| is much closer in spirit to standard statistical simulations of spin models on 
various lattices than the earlier "infinite dimensional" models. It is now well known both from analytical 
work J9, 10 and numerical work jllj that the critical behaviour of spin models such as the Ising model can 



be modified by placing them on appropriate dynamical (connectivity) lattices, where the spins interact 
with the lattice and vice-versa. The methods of || can also be adapted to solve bond percolation on 
dynamical lattices by treating it as the limit of a q state Potts model as q — ► I again leading to 
different critical behaviour from fixed lattice counterparts. In the light of the different behaviour of 
spin models on dynamical lattices and the relation between percolation and damage spreading in the 
finite range Kauffman models a natural question to ask is whether the behaviour of the quenched rule 
automaton is also changed by placing it on a dynamical lattice, which would then correspond to annealing 
the inputs. Such a model lies halfway between the (soluble) model of Derrida and Pomeau with annealed 
rules and inputs and the models with quenched rules and inputs simulated by Derrida and Stauffer, 
and might a priori lie in either universality class. There has been some scepticism of the relevance of 
quenched rule Kauffman automata with nearest neighbour interactions on regular lattices for biological 
applications [ p"3| |, on the grounds that gene interactions are known to have considerable temporal and 
spatial variability. In view of this, our model might even be considered to be as appropriate as the 
original Kauffman infinite range model since we have included both spatial and temporal variation in 
using a dynamical lattice, whilst preserving a characteristic fixed rule for each lattice point or "gene" . 

A second motivation for performing a simulation on a dynamical lattice comes from 2D gravity: in 
113 ] percolation on triangulations arising from simulations of 2D gravity coupled to matter with various 
central charges was used to investigate the geometry of the triangulations. In view of the links between 
damage spreading in the Kauffman model and percolation one might also hope to learn something about 
discretized 2D gravity, at least qualitatively, from simulating the Kauffman model on suitable graphs. 
The geometry of the graphs is certain to be reflected in some of the features of the damage clusters and 
their manner of spreading. There is no back reaction from the automaton spins on the lattice because we 
change the connectivity without reference to the spins so we will, in effect, be investigating the lattices 
of pure 2D gravity with no matter. For the case of pure 2D gravity the appropriate choice of graphs 
would be either <f) 3 graphs or their dual triangulations as in Q . As there is no transition for a quenched 
rule Kauffman automaton on a regular fixed connectivity </> 3 lattice - the honeycomb lattice, we choose 
to simulate 3 graphs to see if the introduction of dynamical connectivity gives a transition. 

In the current work we use the methods employed jll| in our simulations of Potts models coupled 
to 2D gravity where (f> 3 graphs of spherical topology with N — 10000 points were generated. To avoid 
unnecessarily singular graphs, tadpoles and self energy bubbles were not allowed (ie no rings of length 
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one or two in the graphs). The dynamical connectivity was implemented by performing a series of local 
"flip" moves, which are the dual of those employed on triangulations, and can be shown to be ergodic 
Jl5|| . We found that it was a good rule of thumb to carry out NFLIP — N local "flip" moves on the 
lattice in the spin model simulations between spin updates in order to ensure sufficient coupling to 2D 
gravity so we do the same between the automaton time steps. The rules are chosen at random for each 
graph point at the start of each run, using the probabilities p and 1 — p for "up" and "down" rules. As 
the region p > 0.5 is equivalent to p < 0.5 on inversion of p and 1-pwe only simulate up to p = 0.5. The 
fi are then held fixed during the run, whereas the flips ensure that the inputs are shuffled between each 
automaton update, which is performed according to equ.(f ) simultaneously for each spin across the entire 
graph. We look at the evolution of the Hamming distance in equ.(2) for configurations with 4, 100, 400 
and 1000 differing initial spins over 500 automaton timesteps and repeat this process 50 times for each p 
to obtain averaged values and errors. This proved sufficient to obtain reasonable error bars. The initial 
spin configurations were taken as cold (ordered) starts and the initial sites to be damaged were chosen at 
random on each run. We also simulate the automaton on fixed 2D gravity graphs, where we switch off the 
flips in between the automaton updates. In this case we expect to learn something about the geometry 
of a given tfi 3 graph by comparing the results with other fixed graphs, such as the honeycomb graph. 

In Fig.l we have plotted the time evolution of the Hamming distance for p = 0.10 and p — 0.40 
on both fixed and dynamical <fi 3 graphs for a starting damage of 100 to give some idea of the general 
behaviour. One qualitative fact is immediately clear from the graph: there is no transition as p is varied 
for the fixed <p 3 graph, whereas the damage clearly spreads to reach an asymptotic value for p = 0.4 
when the flips are switched on. Thus the fixed 2D gravity graph, a highly disordered <j) 3 graph with a 
complicated internal structure, behaves like its smooth <f) 3 honeycomb counterpart in so far as it displays 
no transition to a chaotic phase at larger p values. The most interesting observation for the dynamical 
graph simulations, on the other hand, is how close they come to the analytical results for the totally 
annealed automaton of Derrida and Pomeau. For a totally annealed Kauffman model with K — 3, which 
is also the value for the dynamical <fi 3 graphs here, one expects a critical value of 

p c = 1/2 - l/(2v / 3) -0.21 (8) 

and an asymptotic value of the damage dfi na i ~ 0.38 at p = 0.5 that is independent of the starting damage. 
In Fig. 2 we plot the asymptotic value of the damage after 500 time steps against p for the various starting 
damages, from which it is clear that the automaton is almost certainly behaving identically to the fully 
annealed model. There is a transition from a frozen to a chaotic phase in the region of 0.21 — 0.22 and 
the final damage is independent of the starting damage, for all but the smallest starting damage, which 
is probably too small to avoid an atypical choice of starting sites. We find dfi na i ~ 0.366(1), which is a 
little lower than the analytical value. The changes in connectivity in the simulation thus appear to be 
sufficient to destroy the correlations between the spin values and allow the application of the annealed 
approximation. 

The results for the asymptotic values of the damage for a fixed 2D gravity graph are shown in Fig. 3, 
clearly demonstrating the absence of a transition in this case. The bond percolation threshold on 2D 
gravity graphs was calculated to be x c ~ 0.78 in so the formula 2p c (l — p c ) = x c of || correctly 
predicts no transition for the quenched Kauffman automaton in this case (as it does for the quenched 
automaton on a honeycomb lattice). It is instructive to compare the results from a similarly sized 
honeycomb lattice in Fig. 4 with those in Fig. 3. The curves for the final damage vs p in Fig. 4 for a lattice 
of 8192 points look, in fact, as if there is a transition. Much larger simulations however, show convincingly 
that a transition is absent on a honeycomb lattice 0, in particular by regarding the behaviour of the 
final damage as the initial damage is varied. From comparing Fig. 3 and Fig. 4 we can therefore deduce 
that it is a lot easier to see the absence of damage spreading on the 2D gravity (j> 3 lattice than on the 
honeycomb (ie regular (j> 3 ) lattice. This suggests that the spreading of a damage cluster encounters more 
obstacles on the 2D gravity graph than the honeycomb. 

Such a conclusion is consistent with the behaviour observed for the critical slowing down of the 
magnetization in cluster algorithms for spin models on dynamical surfaces, which turned out to be 
anomalously large Jll| . It was suggested in jl6| that this was due to the inhibition of cluster growth from 
bottlenecks and "baby-universes" in the graph geometry that only changed slowly due to the local "flip" 



3 



graph update. Such bottle-necks and baby universes would also trap the damage clusters in the quenched 
Kauffman model and account for both the much clearer absence of a transition than on the honeycomb 
lattice and the much smaller final damage on the fixed 2D gravity graphs compared with the honeycomb 
lattice. The relative smallness of the clusters is also an indication of the larger fractal dimension of the 
2D gravity graphs, which appears to be close or equal to 3 O]. 

In summary, we simulate a quenched rule Kauffman automaton on dynamical 3 graphs (as used in 
2D gravity simulations) and find that the annealed results of Derrida and Kauffman are applicable. A 
simulation on a fixed </> 3 graph on the other hand, although there is no transition to a chaotic phase, 
gives some indication of the complicated internal geometry that is present in such graphs. It would 
be interesting to look at the crossover between the effectively annealed behaviour on dynamical graphs 
and that on fixed graphs as the number of flip updates per automaton update was reduced, as well as 
examining the relation between damage spreading and thermodynamic quantities |l8| for more standard 
spin models on dynamical lattices. 
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Figure Captions 

Fig. 1. The damage vs time for p = 0.1,0.4 on both fixed and dynamical 3 graphs with fixed rules. 

Fig. 2. The final damage vs p for various initial damages on a dynamical </> 3 graph with fixed rules. The 
various initial damages are indicated. 

Fig. 3. The final damage vs p for various initial damages on a fixed 4> 3 graph with fixed rules. The 
various initial damages are indicated. 

Fig. 4 The final damage on a similarly sized honeycomb lattice, again with fixed rules. The various 
initial damages are indicated. 
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This figure "figl-l.png" is available in "png" format from: 



http://arXiv.org/ps/hep-lat/9312061vl 
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